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NAVIER-STOKES  CALCULATIONS  WITH  A COUPLED  STRONGLY  IMPLICIT  METHOD 
PART  1:  FINITE-DIFFERENCE  SOLUTIONS 
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Farmingdale,  NY  11735 


tain  a converged  solution.  This  is  a very 
time-consuming  process . 


Stone's  unconditionally  stable,  strongly 
implicit  numerical  method  is  extended  to 
the  2x2  coupled  vorticity  - stream  function 
form  of  the  Navier-Stokes  equations.  The 
solution  algorithm  allows  for  complete 
coupling  of  the  boundary  conditions.  Solu- 
tions for  arbitrary  large  time  steps,  and 
for  cell  Reynolds  numbers  much  greater  than 
two  have  been  obtained.  The  method  con- 
verges quite  rapidly  without  adding  artifi- 
cial viscosity  or  the  necessity  for  under- 
relaxation. This  technique  is  used  here  to 
solve  for  a variety  of  internal  and  exter- 
nal flow  problems.  Moderate  to  large 
Reynolds  numbers  are  considered  for  both 
separated  and  unseparated  flows.  The  pro- 
cedure is  extended  to  higher-order  splines 
in  Part  2 of  this  study. 


In  the  present  paper,  Stone's  strongly 
implicit  method'1  is  reformulated  to  allow 
for  the  solution  of  a 2x2  coupled  system 
of  equations.  The  vorticity-stream  func- 
tion  form  of  the  Navier-Stokes  equations , 
with  second-order  accurate  centered  differ- 
ences, is  solved  by  this  coupled  strongly 
implicit  procedure.  The  solution  algorithm 
also  allows  for  the  coupling  of  the  bound- 
ary conditions  on  the  stream  function  and 
vorticity  in  both  coordinate  directions. 
Solutions  for  arbitrary  large  time  steps 
(At*10°)  and  cell  Reynolds  numbers  much 
greater  than  two  are  possible.  Results 
have  been  obtained  without  adding  artifi- 
cial viscosity  and  without  the  necessity  of 
under relaxation.  Even  with  rather  arbi- 
trary initial  conditions,  the  method  con- 
1 . I ntroduct ion  verges  quite  rapidly  to  the  steady-state 

solution,  when  one  exists. 

A numerical  solution  procedure  for  the 
Navier-Stokes  equations  in  vorticity-stream 
function  form  is  considered.  In  the  past, 
almost  all  explicit  and  implicit  numerical 
systems  that  have  been  developed  for  these 
equations  have  been  solved  iteratively.  It 
has  generally  been  observed  that  for  sta- 
bility purposes,  a nearly  converged  solu- 
tion of  the  Poisson  equation  for  the  stream 
function  is  required  at  each  time  level 
before  incrementing  the  vorticity  equation 
one  time  step.  In  order  to  accelerate  the 
rate  of  convergence  to  the  steady-state, 

SOR  procedures,  modified  ADI  techniques, 
strongly  implicit  methods  and  direct  sol- 
vers have  been  considered  for  the  solution 
of  the  Poisson  equation.  For  the  vorticity 
transport  equation,  explicit  integration 
procedures  are  very  time  consuming  and  the 
temporal  increment  is  restricted  by  the 
CFL  condition.  Implicit  methods,  though 
unconditionally  stable,  usually  exhibit 
spurious  oscillations  and/or  instability 
for  large  time  steps  and/or  large  cell 
Reynolds  numbers.  Moreover,  the  inversion 
of  the  governing  tridiagonal  system  is 
only  assured  if  the  conditions  for  diagonal 
dominance  (Courant  number  less  than  one  and 
cell  Reynolds  number  less  than  two)  are 
satisfied.  The  addition  of  artificial  vis- 
cosity or  the  use  of  upwind  differencing  Finite-difference  solutions  of  partial 

(which  also  contains  large  amounts  of  arti-  differential  equations  of  elliptic  or  para 

ficial  viscosity)  usually  assures  diagonal  bolic  type  require  the  solution  of  alge- 

dominance,  but  lowers  the  over-all  accuracy  braic  systems  of  the  following  type: 
or  adds  excessive  numerical  diffusion. 

Furthermore,  in  many  instances,  it  has  been  A. .W.  . .+B. ,W. .+C.  .W.  . . + 

necessary  to  underrelax  in  the  initial  2 '2  Ip  ^ - r (la) 

stages  of  the  calculation  in  order  to  ob-  °i,j"i-l,j  t,ij*i+l,j  “ij 
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This  procedure  has  been  used  here  to 
solve  for  a variety  of  flow  problems  at 
moderately  large  laminar  Reynolds  numbers, 
generally  R <2000.  The  flow  in  a driven 
cavity,  theeassociated  temperature  field  in 
the  cavity,  the  flow  in  a channel  with  a 
rearward  facing  step,  the  flow  over  a 
finite  length  flat  plate,  and  the  flow  over 
a circular  cylinder  will  be  presented. 

In  the  next  section,  the  solution  proce- 
dure is  described.  One  of  the  important 
limitations  of  the  present  technique  is  the 
comparatively  large  storage  requirement  for 
the  inversion  algorithm.  However,  in  the 
absence  of  any  direct  solvers  for  the  non- 
linear coupled  equations  considered  here, 
the  ability  to  obtain  converged  Navier- 
Stokes  solutions  at  moderately  large 
Reynolds  numbers  in  modest  computer  times 
is  particularly  attractive.  In  a second 
part  of  this  study,  higher-order  spline 
collocation  procedures  are  applied  with 
this  coupled  algorithm.  This  significantly 
reduces  the  grid  and  therefore  the  storage 
requirements . 


In  matrix  form  this  can  be  written  aa 


where  L is  a NxN  matrix  of  the  following 
form: 


leads  to  the  familiar  point  Jacobi-itera- 
tive procedure.  The  Gauss-Siedel*  method 
is  recovered  by  chosing  P such  that  A is, 
the  lower  triangular  part  of  the  matrix  L. 

Although  many  other  variants  are  possible, 
the  primary  contribution  of  Stone  has  been 
to  devise  the  factorization  in  such  a way 
that  a certain  degree  of  implicitness  is 
associated  with  each  coordinate  direction 
and  such  that  every  element  in  PW^  is 
small;  in  particular,  the  elements  are  of 
order  h , where  h is  the  mesh  width.  How- 
ever, Saylor”  infers  that  a first-order 
factorization  can  be  more  useful  than  those 
that  lead  to  second-order  correction  terms. 
In  any  event,  PW"  should  tend  to  zero  as  h 
vanishes.  An  undesirable  feature  of  Stone's 
factorization  is  that  the  matrix  P changes 
at  each  step  so  that  two  successive  itera- 
tions are  in  a sense  uncorrelated.  Al- 
though, the  ideal  factorization  may  depend 
upon  the  particular  problem  being  con- 
sidered or  upon  previous  experience,  the 
algorithm  presented  herein  is  sufficiently 
general  for  all  types  of  factorization. 

The  final  two-pass  algorithm  can  be  written 
as: 

(L+P)  Wn+1  = AW11*1  = G+Pw"  (3) 


The  solution  of  these  equations  have  been 
obtained  by  a variety  of  techniques,  in- 
cluding direct  solution  by  Gaussian  elimin- 
ation or  iterative  methods  such  as  SOR, 

ADI,  etc.  Stone1  developed  a strongly 
iiuplicit  method,  which  is  more  rapidly  con- 
vergent and  is  based  on  a quasi-LU  decom- 
position of  the  matrix  L. 

Direct  solvers  such  as  those  due  to  ^une- 
man,  Sweet  and  Schwarztrauber,  Bank, 
etc.  are  only  applicable  to  a special 
class  of  the  difference  equations  (1)  and 
associated  boundary  conditions.  These  di- 
rect methods  are  not  presently  useful  for 
the  general  finite-difference  form  of  each 
of  the  Navier-Stokes  equations  and  certain- 
ly not  for  the  coupled  vorticity-stream 
function  system.  With  appropriate  boundary 
conditions  the  Poisson  equation  is  amenable 
to  direct  methods;  however,  some  iterative 
numerical  procedure  is  required  for  the  un- 
coupled vorticity  equation.  Other  variants 
of  Gaussian  elimination  are  extremely  inef- 
ficient and  time  consuming  and  even  sus- 
ceptible to  a large  accumulation  of  round- 
off error.  SOR  and  ADI  techniques  converge 
rather  slowly  and  for  certain  types  of 
differencing  (involving  9 points)  SOR  may 
not  converge  at  all. 

Stone's  strongly  implicit  iterative  tech- 
nique falls  under  the  general  category  of 
“factorization  methods".  The  underlying 
idea  of  factorization  is  to  replace  the 
sparse  matrix  L by  a modified  form  (L+P) 
such  that  the  resulting  matrix  can  be  de- 
composed into  upper  (0)  and  lower  (L)  tri- 
angular sparse  matrices.  This  leads  to  the 
following  general  iterative  procedure  for 
the  system^  (1) : 

(L+P),Wn+1  - G+PW"  , (3a) 


LV  = G+PW"  , UWn+1  = V (4a' b) 

In  Stone's  factorization  procedure  P is 
prescribed  such  that  L and  U have  only 
three  non- zero  diagonals.  This  leads  to  a 
solution  of  the  following  form: 

«if  = + Eij  wi,S+i + Fij  wi+f  j.(4c) 

This  procedure  has  the  distinct  advantage 
of  being  implicit  in  both  the  i and  j di- 
rections, as  well  as  coupling  all  the 
boundary  conditions.  We  have  found  that 
this  technique  generally  converges  more 
rapidly  than  do  many  of  the  more  familiar, 
less  implicit,  iterative  methods  previous- 
ly mentioned.  Moreover,  it  is  the  present 
authors’  opinion  that  the  lack  of  implicit- 
ness and  coupling  may  explain  some  of  the 
difficulties  that  others  have  encountered 
in  their  solutions  of  high  Reynolds  number 
Navier-Stokes  flows,  in  particular,  those 
concerning  the  CPL  limitation,  the  need 
for  artificial  viscosity  and  the  need  for 
underrelaxation.  In  view  of  these  observa- 
tions it  would  appear  that  a direct  solver 
would  be  siost  suitable  for  the  solution  of 
the  algebraic  system  of  difference  equa- 
tions (1)  arising  from  the  Navier-Stokes 
equations.  Unfortuantely,  as  noted  pre- 


then 


where  the  superscript  n denotes  the  itera- 
tion number.  The  rate  of  convergence  of 
this  iteration  scheme  depends  upon  the 
particular  ehoicA  of  the  matrix  P.  The 
two  essential  requirements  on  the  matrix  P 
are  as  follows:  i)  the  elements  of  P 
should  be  small  in  magnitude  so  that  the 
explicit  perturbation  is  small  and  ii)  the 
resulting  matrix  A should  be  decomposable 
into  sparse  L,U  factors.  If  the  L and  U 


viously,  currently  available  efficient 
direct  solvers  are  not  applicable  to 
coupled  2x2  non-linear  systems,  while 
others  are  inefficient,  time  consuming  and 
suffer  from  instability  due  to  round-off 
accumulation.  The  strongly  implicit  itera- 
tive procedure  can  provide  the  necessary 
coupling  of  the  dependent  variables  and 
boundary  conditions,  and  add  a significant 
degree  of  efficiency  and  speed  of  con- 
vergence. On  the  other  hand,  computer 
storage  requirements  are  increased.  This 
problem  can  be  minimized  by  reducing  the 
required  number  of  grid  points.  With  the 
higher-order  methods  to  be  presented  in 
Part  2 this  becomes  possible. 

In  the  next  section,  the  factorization 
procedure  will  be  developed  for  a coupled 
2x2  system,  e.g.,  the  vorticity-stream 
function  form  of  the  two-dimensional  Navier 
Stokes  equations. 

3.  Coupled  2x2  Algorithm 

Jacobs^ hae  extended  Stone's  algorithm  for 
the  second-order  accurate  nine-point  cen- 
tral difference  form  of  the  biharmonic  op- 
erator. He  applied  the  algorithm  to  the 
stream  function  form  of  the  Navier-Stokes 
equations.  The  algorithm  was  used  to  in- 
vestigate the  flow  in  a driven  cavity. 
Jacobs could  not  obtain  converged  solutions 
for  large  Reynolds  numbers,  i.e.,  R >400. 
The  failure  of  the  algorithm  could  Be 
associated  with  the  iterative  procedure  for 
implementing  the  boundary  conditions.  Our 
own  experience  with  the  strongly  implicit 
procedure  indicates  that  the  use  of  itera- 
tive boundary  conditions  considerably  re- 
duces the  applicability  (convergence  rate, 
stability)  of  the  algorithm.  Therefore,  a 
2x2  coupled  strongly  implicit  algorithm 
has  been  developed  for  the  stream-function/ 
vorticity  form  of  the  Navier-Stokes  equa- 
tion. The  finite-difference  form  of  these 
equations  are  as  follows: 


The  iterative  solution  is  given  as 
IL+P]  Wn+l  « G + pw" 


As  described  in  Section  2,  the  following 
two-pass  algorithm  for  equations  (5a, b)  is 
obtained: 


See  the  Appendix  for  the  recurrence  rela 
tionships . 


Note  that  the  evaluation  of  u>.  . depends 
not  only  upon  the  other  values  of  id,  but 
also  upon  the  new  values  of  t|>.  The  coup- 
ling of  id  and  ip  at  the  boundaries  is  ex- 
act. In  a single  iteration,  the  calcula- 
tion along  arrays  close  to  the  boundaries 
is  very  close  to  an  exact  solution.  When 
the  boundary  conditions  (no-slip  condi- 
tions) are  satisfied  iteratively,  this  fea 
ture  of  the  algorithm  is  lost  and  the  solu- 
tion converges  more  slowly.  In  these  it- 
erative procedures,  the  calculation  may 
even  become  unstable.  It  is  for  these  rea- 
sons that  underrelaxation  is  sometimes  re- 
quired. 

The  uncoupled  algorithm,  equation  (4) , 
requires  three  auxiliary  functions;  the 
present  2x2  algorithm  requires  10  auxili- 
ary functions.  This  represents  an  in- 
crease in  storage  requirements  by  a factor 
of  more  than  three.  Unfortunately,  this 
is  one  of  the  serious  drawbacks  of  the 
present  method.  However  with  the  (2x2) 
coupled  method,  steady-state  calculations, 
without  CFL  limitations,  artificial  vis- 
cosity, or  underrelaxation,  can  be  per- 
formed in  relatively  few  iterations.  With 
a smaller  computational  field,  i.e.,  a 
coarse  mesh,  the  storage  requirements  *re 
reduced.  Therefore,  higher-order  tech- 
niques for  obtaining  steady-state  solu- 


Szcflu 

Section 


Direct  solution  of  these  equations  by  var- 
iants of  Gaussian  elimination  is  very  in- 
efficient and  expensive.  The  strongly  im- 
plicit procedure,  though  iterative  in 
nature,  provides  the  appropriate  coupling 
and  implicitness  required  for  convergence 
The  factorization  procedure  given  in  the 
preceding  section  is  quite  general  and  is 
applicable  to  equations  (5a, b)  as  well. 
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1 
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tions  with  the  present  Method  become  very 
desirable.  Higher-order  spline  4 solutions 
of  the  Navier-Stokes  equations  with  the 
2x2  coupled  algorithm  are  presented  in 
Part  2 of. the  present  paper.  In  earlier 
studies  ' of  the  spline  4 procedure  for 
the  Navier-Stokes  equations,  a reduction  of 
a factor  of  sixteen  in  the  number  of  grid 
points  has  been  realized  in  order  to  obtain 
accuracy  comparable  with  the  second-order 
central  difference  solution.  This  reduces 
the  over-all  storage,  required  for  execu- 
tion of  the  algorithm,  by  a factor  of  160. 


Governing  Equat ions 


The  vorticity-stream  function  form  of  the 
two-dimensional  Navier-Stokes  equations  are 
considered  for  a number  of  internal  and  ex- 
ternal flows.  The  calculations  are  carried 
out  in  a transformed  (x,y)  coordinate 


system  using  conformal  mappings  to  simplify 
the  computational  domain. 


j -j^r  ♦ <u*>x  ♦ (vw)^  « g—  V2w  (7a) 


where  R = Ua/v;  U,a  are  reference  velocity 
and  length;  v is  the  kinematic  viscosity; 

Pr  is  the  Prandtl  number;  u>  and  4 are  the 
vorticity  and  stream  function;  u and  v are 
the  velocities  along  the  x and  y coordin- 
ate directions  in  the  transformed  plane;  T 
is  the  temperature  and  J is  the  Jacobian  of 
the  transformation. 


'he  mapping 


is  considered,  where  J = and 

Z = x\+  iy,  C * x + iy.  s (x,y)  are  the 
physical  coordinates.  It  should  be  noted 
that  tne  velocities  (u,v)  in  the  trans- 
formed_blane  are  not  the  'physical  veloci- 
ties (u/y)  in  the  (x,y)  directions.  The 
mapping  defines  the  relationship  between 
(u,v)  and\(u,v). 


Although \ the  equations  (8)  are  written 
for  unsteady  flow,  only  the  steady-state 
is  considered  in  the  present  study,  there- 
fore very  large  temporal  increments 
(At~10°) are  prescribed.  As  noted  previous' 
ly,  this  is  pcmsible  with  the  coupled 
strongly  implicit  procedure.  The  appropri' 
ate  boundary  conditions  on  the  solid  sur- 
faces are  the  nq-slip  conditions  and  the 
Dirichlet  condition  for  the  stream  func- 
tion. As  described  by  Roache,  the  calcu- 
lations can  be  veW  sensitive  to  the  in- 
flow and  outflow  boundary  conditions;  in 
the  present  paper , 'derivative  boundary  con 
di tions  have  been  utilized  at  these  bound- 
aries. These  will  tfe  detailed  for  the 
problems  to  be  considered  here. 


5.  Pinite-Pifference  Discretization 

The  following  difference  expressions  are 
prescribed  for  the  derivatives  $x,  $xx: 

♦x  * * ki-xi-xi-l  ' °i-ki+lAi 

2(*i+l-(1+'°l,*i+0l*i-l)  <8) 

_ p 

xx  a^l+o^kf 

where  4 is  any  of  the  physical  variables. 

A non-uniform  grid  is  prescribed  with 
a.  = 1+0 (k. ) in  order  to  maintain  second- 
order  accuracy  of  the  difference  formulas. 
If  u and  v,  in  the  vorticity  transport 
equation,  are  treated  explicitly,  the 
governing  equations  are  implicitly  coupled 
only  through  the  boundary  conditions.  In 
order  to  strengthen  the  coupling,  the 
non-linear  convective  terms  will  be  quasi: 
linearized.  The  K-R  differencing  scheme10 
is  applied  for  the  vorticity  transport 
equations.  This  differencing  technique  is 
in  effect  a splitting  of  the  convective 
terms  into  an  implicit  upwind  term  and  an 
explicit  effective  diffusive  corrector. 
Physical  diffusion  terms  are  treated  im- 
plicitly. The  converged  steady-state 
solution  is  second-order  accurate.  This 
deferred  correction  procedure  although  ap- 
parently obvious,  has  not  generally  been 
used  and  has  been  found  to  be  of  major  im- 
portance not  only  in  the  present  context 
but  in  the  application  of  higher-order 
numerical  methods.  The  extension  of  this 
scheme  to  fourth  and  sixth-order  methods, 
and  the  associated  stability  analysis  has 
been  presented  in  reference  (9) . With  the 
K-R  scheme,  the  convective  terms  are  writ- 
ten as: 

(uto)i  - ,-<uu>).  .n+1 

(U“>X  - °X  1 ^ ^ + 

(uid)  . . - (uu>)  . . .n+1 

♦ ‘1-»*x>  f — LfJ)  + 


+ {1-M  (1+i-)  )Dn  , 

A U ^ 


where  D = 


(UO>)  i + j-  (ltQ1)  (Ufa))  (uin)  ■ 


and  ux  * 0 ifu^ ^>0  ; ux  = 1 ifui^<0 

Similar  expressions  are  used  for  (vw)  . 

The  terms  (uo>)  . . and  (vw.J  are  quasi-y 
linearized  to  provide  the^necessary  5-  « , 

point  coupling  between  4 and  u during  the 
iterative  procedure.  It  should  be  noted 
that  there  are  various  other  methods  for 
achieving  this  coupling,  e.g.,  by  adding 
and  subtracting  the  upwind  differencing 
terms  at  the  n+1  and  n time  levels,  re- 


for  (7a,b)  is  independent  of  the  coupling 
procedure  and  the  rates  of  convergence  were 
not  noticeably  different  for  the  problems 
considered  here.  The  ideal  splitting  for 
the  convective  terms  may  depend  upon  the 
method  of  solution  of  the  algebraic  system 
and  should  provide  the  maximum  rate  of  con- 
vergence. The  energy  equation  (7c)  is 
solved  uncoupled  from  the  (ij>,(i>)  system  once 
the  u, v values  have  been  obtained. 

6.  Boundary  Conditions 

The  boundary  conditions  on  i p and  u for 
the  strongly  implicit  procedure  are  satis- 
fied in  the  same  manner  as  for  the  simple 
one-dimensional  LU  decomposition  or  as  it- 
is  sometimes  termed  the  Thomas  algorithm. 
The  distinct  advantage  of  the  present 
coupled  procedure  is  that  all  the  bounda- 
ries interact  implicitly  and  simultaneous- 
ly. For  example,  if  we  consider  the  no- 
slip condition  at  a solid  surface,  the 
usual  first-order  boundary  condition  for 
(<p,u)  is 


tions.  Underrelaxation  was  reauired  in  the 


initial  stages  of  the  calcuation.  A con- 
verged solution  for  the  stream  function  is 
obtained  and  then  the  vorticity  is  updated 
in  the  time  iteration.  Many  iterative 
studies  have  followed;  more  recently, 
Buneman's  direct  solver  for  ij/  and  the  hop- 
scotch or  ADI  methods  for  the  conservation 
form  of  the  vorticity  transport  equation 
were  employed  successfully  by  Smith,  J He 
has  been  able  to  obtain  central  difference 
solutions  for  R <5000.  This  procedure  re- 
quires a small  St  (though  not  limited  by 
the  CFL-condition)  for  stability  and  con- 
sequently takes  a longer  time  for  conver- 
gence to  the  steady  state.  Nallasamy  et 
al.  applied  upwind  differencing  to  ob- 
tain solutions  for  R <50,000.  This  proce- 
dure however  suf ferseTrom  large  numerical 
viscosity  and  therefore  the  effective  R 
is  in  fact  very  low.  The  application  of 
the  strongly  implicit  method  for  the  bi- 
harmonic stream  function  form  of  the 
Navier-Stokes  equations  was  considered  by 
Jacobs.  As  previously  noted,  these  cal- 
culations diverged  for  Re>400. 

In  the  present  paper,  solutions  for, 

R <3000  have  been  obtained  with  At=10  , a 
17x17  grid,  and  with  the  strongly  implicit 
coupled  procedure  outlined  previously. 
Arbitrary  initial  conditions  were  assumed. 
The  results  obtained  here  are  identical  to 
those  found  previously  by  direct  solver, 
SOR  or  ADI  methods,  when  such  solutions 
were  possible.  The  velocity  distribution 
through  the  vortex  center  for  R =1000  and 
2000  is  shown  in  Fig.  2,  the  values  of 


In  an  uncoupled  procedure,  the  wall  vorti- 
city is  updated  after  <p  has  been  evaluated 
everywhere.  In  the  present  coupled  proce- 
dure, equation  (10)  is  satisfied  exactly 
in  each  iteration.  In  order  to  fully 
utilize  the  positive  features  of  the  al- 
gorithm, it  is  highly  desirable  that  all 
the  boundary  conditions  are  coupled  im- 
plicitly. The  second-order  no-slip  con- 
dition (11)  is  used  for  all  of  the  present 
calculations 


,*<l+0-)l|l.  5+0,  <l> 


The  bracketed  portion  is  treated  implicit- 
ly. The  correction  term  enters  explicitly 
The  deferred  approach  for  the  higher-order 
correction  to  the  boundary  condition  has 
also  been  applied  for  the  central  differ- 
ence corrector  to  upwind  differencing  in 
what  is  termed  the  K-R  procedure10  and  for 
higher-order  splines  in  Part  2 of  the  pres' 
ent  analysis,  see  reference  (9). 


7.1  Cavity  Flow 


The  flow  in  a driven  cavity,  see  Fig.  1, 
has  been  investigated  by  many  authors 
using  a variety  of  iterative  procedures:, 
e.g.,  see  references  (11-13).  Burggraf1* 
in  one  of  the  first  analyses  obtained  solu' 
tions  for  R < 400  with  the  SOR  method  and 
the  steady  Kate  form  of  the  (4i-w)  equa- 


7.3  Flow  in  a Channel  with  a Backward  Facin' 


maximum  stream  function,  position  of  the 
vortex  center  and  the  value  of  the  vorti- 
city  are  also  depicted.  In  addition. 

Smith's  central  difference  solution,  but 
with  first-order  accurate  boundary  condi- 
tions are  also  shown.  It  should  be  noted 
that  the  accuracy  of  the  wall  vorticity  has 
a marked  influence  on  the  solution.  The 
present  method  converges  in  50  to  100  it- 
erations as  R changes  from  100  to  3000. 

For  the  same  grid.  Smith  required  about  120 
and  693  iterations  for  R =100  and  1000,  re- 
spectively. The  number  of  iterations  re- 
quired by  Smith  increases  to  3469  for 
R =5000.  Though  the  application  of  the 
strongly  implicit  coupled  procedure  repre- 
sents a considerable  improvement  over  the 
time  dependent  methods,  the  rate  of  con- 
vergence associated  with  the  present  zeroth- 
order  factorization  procedure  can  still  be 
improved.  Various  methods  for  increasing 
the  rate  of  convergence  will  be  examined  in 
a future  study. 

Summarizing  our  conclusions  for  the  cavity 
problem,  the  time  dependent  methods  require 
large  numbers  of  iterations  even  with  a 
direct  solver  for  the  stream  function.  Up- 
wind differencing  results  must  be  inter- 
preted with  caution  due  to  the  presence  of 
large  amounts  of  artificial  viscosity.  In 
certain  cases,  as  will  be  shown  for  the 
flow  over  a circular  cylinder,  upwind  dif- 
ferencing can  lead  to  completely  incorrect 
conclusions.  The  present  coupled  method, 
though  reasonably  fast  does  require  con- 
siderably more  computer  storage  than  the 
other  techniques  for  large  numbers  of  grid 
points.  This  problem  decreases  with  the 
higher-order  spline  methods  to  be  presented 
in  Part  2. 

7.2  Heat  Transfer  in  a Driven  Cavit 


Ml  NEXT  TRANSFER  ON  SIOE 
OF  CANITY 


This  problem  was  considered  as  one  of  the 
test  cases  at  the  International  Conference 
on  Numerical  Methods  in  Laminar  and  Turbu- 
lent Flows  held  at  University  College  of 
Swansea.  The  (u,v)  velocity  distribution 
was  prescribed.  The  origin  of  these  pro- 
files is  unknown. 


In  order  to  satisfy  the  adiabatic  condition 
to  second-order,  the  Tyy  term  at  the  wall 
was  evaluated  by  incorporating  the  govern- 
ing equation  (7c) . Once  again,  all  the 
boundary  conditions  are  implicitly  coupled 
into  Stone's  solution  algorithm.  The  en- 
ergy equation  is  differenced  in  non-con- 
servation form.  The  solution  for 
P =R  .p  =50  is  obtained  on  a 15x15  uniform 
gfid?  fhe  solutions  for  the  heat  transfer 
at  the  side  walls  is  shown  in  Fig.  3;  the 
temperature  distribution  on  the  upper  and 
lower  walls,  as  well  as  the  mid  plane,  is 
shown  in  Fig.  4.  It  should  be  noted  that 


ns.  4 TEMPERATURE  PROFEES  IN  CANITY. 


the  prescribed  velocity  distribution  does 
not  exactly  satisfy  the  continuity  equation 
as  applied  herein.  Therefore,  the  present 
results  are  only  as  good  as  the  accuracy  of 
these  velocity  conditions . 


The  Navier-Stokes  equations  in  primitive 
variables  have  been  solved  by  Taylor  and 
Ndefo16  for  the  step  geometry.  Solutions 
for  R =25,  100  were  obtained, using  explicit 
techniques . Recently  Roache  has  investi- 
gated the  equations  for  this  problem. 

Steady  state  equations  were  considered; 
however,  underrelaxation  was  required  to 
stabilize  the  solution  for  the  vorticity 
transport  equation.  Also,  this  method  of 
solution  is  restricted  to  large  grid  aspect 
ratio. 


In  the  present  study,  the  2x2  coupled  al' 
gorithm  has  been  used  to  solve  the  4i-o) 
equations  for  the  step  channel.  Steady 
state  (At=10°)  solutions  for  R =1,  100  and 
1000  have  been  obtained.  The  reference 


The  thermal  boundary  conditions  are  shown 
in  Fig.  1.  The  adiabatic  wall  boundary 
condition  is  satisfied  by  a Taylor  series 
expansion  of  the  form: 


Ti,2”Ti,l+h(Ty)i,l+5-(Tyy)i,l+0(h3)  ’ (12) 
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length  a in  R is  the  channel  width  k.  U 
is  the  mean  inflow  velocity.  The  step  ge- 
ometry, shown  in  Fig.  5 is  first  trans- 
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formed  into  a straight  channel  by  using  the 
following  conformal  transformation: 


„_k v-1 ,2w-c-l,  k ,.-1 , (c+])w-2c  ,,,v 

Z-^cos  h ( c-T  )-- -cos  h ( 

IT  /C 


log  w. 


He 


1/2 


where  H is  the  step  height  and  k the  width 
of  the  channel.  A 48x21  variable  grid, 
see  Fig.  6,  is  employed  in  the  c plane. 

The  mesh  aspect  ratio  varies  from  a very 
small  number  (0.03)  to  a large  number 
(5810) . The  convergence  of  the  solution 
procedure  was  not  sensitive  to  these  large 
changes  in  the  grid  aspect  ratio.  At  the 
inflow,  fully  developed  conditions  are  as- 
sumed. These  conditions  are  assumed  in 
order  to  avoid  the  effect  of  the  entrance 
boundary  layers  on  the  recirculating  flow 
region  behind  the  step.  At  the  outflow, 
the  boundary  layer  equations  are  used  im- 
plicitly as  a part  of  the  coupled  al- 
gorithm. Fully  developed  flow  conditions, 
which  change  discontinuous ly  at  the  step 
were  assumed  for  the  initial  guess.  It 
should  be  pointed  out  that  under relaxation 
or  artificial  viscosity  were  not  required 
for  convergence.  For  the  various  Reynolds 
numbers  considered  here,  convergence  was 
achieved  in  from  40  to  80  iterations;  this 
compares  with  2000  required  by  Taylor  and 
Ndejo. 


In  the  present  calculations,  a Reynolds 
number  scaling  along  the  channel  was  not 
imposed,  to  provide  optimal  resolution  in 
the  recirculation  bubble.  The  region  of 
separated  flow  increases  with  Reynolds 
number  and  this  is  evident  from  the  pres- 
ent results.  The  distribution  of  vorti- 
city  along  the  upper  and  lower  walls  is 
shown  in  Figs.  7,8  for  R *1,  100,  1000. 

The  vorticity  near  the  step  is  shown  on  an 
expanded  scale  in  Fig.  8.  Vector  plots  of 


FIG.  6 GRD  DISTRIBUTION  FOR  CHANNEL 
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the  streamlines  are  shown  in  Fig.  9.  It 
can  be  seen  that  for  R *1,  the  flow  separ- 
ates below  the  corner.  As  the  Reynolds 
number  increases,  the  separation  point 
moves  upward  toward  the  corner.  Since  the 
transformation  (13)  is  singular  at  the 
corner,  a grid  point  was  not  located  at 
the  juncture  itself;  however,  from  the 
calcuated  results,  the  indications  are 
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fig. 9b  flow  vector  ougram:  channel-  R,»  100 

that  for  larger  Reynolds  numbers  the  separ- 
ation point  moves  even  closer  to  the  cor- 
ner. Although  the  inflow  profiles  are  in- 
dependent of  the  Reynolds  number,  the  flow 
approaching  the  corner  acquires  the  charac- 
teristics typical  of  a high  R boundary 
layer.  The  degree  of  upstreafi  influence  is 
an  inverse  function  of  the  Reynolds  number. 

7.4  Flow  Over  a Flat  Plate 

Both  semi-infinite  and  finite  flat  pla§e 
flows  for  unit  Reynolds  numbers  up  to  10 
have  been  obtained  with  the  2x2  coupled  al- 
gorithm. Step  initial  condition  were  as- 
sumed. These  were  applied  at  the  leading 
edge  for  the  semi-infinite  plate  so  that  no 
upstream  influence  was  allowed.  This 
to  velocity  overshoots  near  the  inflow  which 
slowly  decay  downstream.  On  the  same  grid 


FIG  9c  FLOW  VECTOR  DIAGRAM'.  CHANNEL- R§*  IOOO 


a solution  of  the  ili-w  form  of  the  boundary- 
layer  equations  was  also  obtained.  The  two 
results  compared  very  closely.  Therefore 
the  relative  accuracy  and  consistency  of 
the  central-difference  schemes  for  the 
Navier-Stokes  equation  and  very  large  R 
has  been  demonstrated.  If  a boundary- lSyer 
solution  is  used  at  the  inflow  instead  of 
the  step  conditions,  the  overshoots  in 
velocity  disappear.  Calculations  for  the 
finite  flat  plate  with  uniform  flow  far  up- 
stream of  the  leading  edge  were  also  con- 
vergent and  acceptable. 

7 . 5 Flow  Over  a Circular  Cylinder  with  a 


described  by  Dennis  and  Chang  for  R =100 
The  question  of  a steady-state  limit  ap- 
pears to  remain  unresolved. 


In  the  present  study  the  coupled  2x2  al- 
gorithm is  used  to  reexamine  the  existence 
of  steady  flow  solutions.  The  cylinder 
with  a splitter  plate  is  mapped  into  a 
semi-infinite  flat  plate  with  the  Joukowski 
transformation,  see  Fig.  10. 


The  numerical  evaluation  of  the  flow  over 
a circular  cylinder  has  been  the  subject  of 
investigation  of  many  authors. -.Thom, 
Kawaguti,..  Takami  and  Keller,  and  Dennis 
and  Chang  have  considered  the  steady 
state  equations  for  The  resulting 

second-order  central  difference  equations 
were  solved  by  SOR  or  line-relaxation  tech- 
niques,. Except  for  the  work  of  Dennis  and 
Chang,  who  obtained  solutions  for  R =100, 
these  steady-state  calculations  have  Been 
limited  to  R <60.  The  reference  length  in 
Re  is  the  cylinder  diameter. 

The  preponderance  of  studies  have  used 
time -dependent  methods  to  evaluate  the 
transient  and  asymptotic  steady  state  be- 
havior. For  large  R flows,  the  solution 
for  the  wake  region  Is  apparently  oscilla- 
tory and  the  question  of  the  existence  of  a 
steady  state  solution  remains  open.  Almost 
all  time-dependent  calculations  suggest 
that  these  wake  flow  oscillations  appear 
for  R >100.  Recently,  Lin  et.  al.  rein- 
vestifated  the  cylinder  flow  with  the 
strongly  implicit  method.  The  i|>-u>  equa- 
tions were  solved  iteratively  using  Stone's 
algorithm  described  herein.  Oscillatory 
wake  flow  was  observed  for  R >80.  This 
contrasts  with  the  "steady-state'*  solution 


PM.  10  CYLINDER  GEOMETRY 

Only  the  upper  half  of  the  flow  was  con- 
sidered. A 33x17  variable  grid  generated 
from  an  analytical  transformation  was 
prescribed.  This  grid  is  very  coarse,  but 
adequte  for  the  present  purposes.  The 
mesh  width  at  the  surface  is  (Ay) 2-25. 
Fifteen  of  the  33  points  in  the  axial  di- 
rection are  on  the  surface  of  the  cylinder 
see  Fig.  10.  Derivative  boundary  condi- 
tions 


<*>  * °,  <Py  = 1 as  y-*® 

“xx  = 0 = ^xx  as  <15) 

were  employed  for  the  upper  and  outflow 
boundaries  respectively.  For  R <70  con- 
verged solutions  were  obtained  for  the 
steady  state  equations.  Solutions  for  the 
wall  vorticity  for  R =50,  60,  70  are  shown 
in  Fig.  11.  For  R =70  the  separation  oc- 
curs at  6=61.8°.  for  R =80  and  100,  a 
steady  flow  in  the  wakeeregion  was  not  ob- 
tained. For  these  values  of  Rg,  a con- 


mu  voffnarr  along  surmcc  or  cyunogr 

verged  steady  state  solution  for  the  up- 
wind difference  equations  was  obtained. 

The  second-order  correction,  using  the  un- 
conditionally stable  K-R  scheme  did  not 
converge  for  At=10.  Solutions  were  found 
for  smaller  values  of  At.  The  solutions 
for  the  wall  vorticity  for  R =100  and 
A t=. 1,1, 10  are  shown  in  Fig.  12.  The  cal- 


culations for  this  problem  were  carried 
out  on  PDP  11/60  in  an  interactive  mode. 
The  solution  was  monitored  at  each  time 
step.  Even  larger  values  of  At (=100) 
were  used  without  any  difficulty.  The 
solution  however  never  acquired  a steady- 
state  value.  As  seen  from  Fig.  12,  al- 
though the  vorticity  in  the  wake  chances 
considerably,  the  effect  on  the  solution 


fore  and  aft  of  the  cylinder  recirculation 
region  is  negligibly  small.  Moreover,  the 
oscillatory  behavior  appears  to  be  closely 
related  to  movement  of  the  point  of  reat- 
tachment in  the  wake.  This  non-stationary 
behavior  is  present  even  with  the  splitter 
plate.  Although  the  splitter  plate  sup- 
presses the  antisymmetric  Karman  vortex 
street,  it  does  not  seem  to  eliminate  the 
vortex  shedding  completely  and  therefore 
the  solution  remains  only  quasi-steady. 

Significantly,  steady  state  upwind  dif- 
ference solution,  even  for  R =1000  and 
At=10  , were  obtained.  The  Converged  solu- 
tion for  wall  vorticity  is  presented  in 
Fig.  13.  It  should  be  noted  that  the  pres- 
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ence  of  large  amounts  of  artificial  viscos- 
ity for  the  very  crude  grid  considered  here 
reduces  the  size  of  the  separated  flow 
region  suggesting  a much  smaller  effective 
R . This  allows  for  a steady  state  solu- 
tion, for  a flow  which  in  fact  is  oscilla- 
tory. In  the  light  of  this  result  and 
others  found  for  the  cavity  flow,  upwind 
difference  results  must  be  viewed  with 
great  caution,  and  in  the  opinion  of  the 
present  authors  may  not  be  meaningful  un- 
less the  mesh  size  is  rather  small. 

8 . Summary 

The  incompressible  vorticity-stream  func- 
tion system  of  the  Navier-Stokes  equations 
has  been  solved  numerically  for  a variety 
of  flow  configurations:  (i)  driven  cavity, 
(ii)  thermal  cavity,  (iii)  channel  with  a 
step,  (iv)  flat  plate,  (v)  circular  cylin- 
der with  a splitter  plate.  The  important 
features  of  the  present  analysis  are  as 
follows:  1)  The  coupled  (w-i (i)  equations, 

with  coupled  boundary  conditions,  are  con- 
sidered. 2)  A strongly  implicit  procedure 
is  applied.  3)  The  combination  of  1)  and 
2)  allows  for  solutions  at  relatively  large 
laminar  Reynolds  numbers  (R  <3,000);  flows 
with  separation  regions  areeincluded. 

4)  The  time  step  can  be  taken  arbitrarily 
large  so  that  in  effect  this  is  a steady- 
state  solver.  Convergence  is  achieved  in 
30  to  100  iterations.  5)  Large  spatial 
grids  are  allowable  so  that  cell  Reynolds 


numbers  are  considerably  in  excess  of  two. 
6)  No  artificial  viscosity  is  required  to 
obtain  second-order  accurate  central-dif- 
ference solutions.  7)  For  the  cylinder 
geometry,  steady  solutions  were  found  only 
for  R <80.  For  larger  Reynolds  numbers, 
transient  calculations  for  smaller  time 
steps  depict  an  oscillatory  behavior,  that 
is  essentially  confined  to  separation 
region.  8)  It  has  been  demonstrated  that 
strongly  implicit  solutions  are  possible 
for  all  R examined  with  upwind  differ- 
encing. This  is  a result  of  the  large 
numerical  viscosity,  which  lowers  the  ef- 
fective Rfi  significantly. 

References 


1.  Stone,  H.L.:  "Iterative  Solution  of 
Implicit  Approximations  of  Multidimension- 
al Partial  Equations".  SIAM  J.  Numer. 
Anal.,  Vol.  5,  pp.  530-558  (1968). 

2.  Roache,  P.J.:  "Computational  Fluid  Dy- 
namics".  Hermosa  Publishers  (1972). 

3.  Buneman,  O. : "A  Compact  Non-iterative 
Poisson  Solver".  SUIPR  Report  294  Inst, 
for  Plasma  Research,  Stanford  Univ.  Calif. 
(1969)  . 

4.  Schwarztrauber,  P.N.  and  Sweet,  R.A.: 
"The  Direct  Solution  of  the  Discrete 
Poisson  Equation  on  a Disc".  SIAM  J. 

Numer.  Anal.,  Vol.  10,  pp.  900-907  (1973). 

5.  Bank,  R.E.:  "Marching  Algorithms  for 
Elliptic  Boundary  Value  Problems.  II: 

The  Variable  Coefficient  Case:.  SIAM  J. 
Numer.  Anal.,  Vol.  !5,  pp.  950-970  (1977). 

6.  Saylor,  P.:  "Second  Order  Strongly 

Implicit  Symmetric  Factorization  Methods 
for  the  Solution  of  Elliptic  Difference 
Equations:.  SIAM  J.  Numer.  Anal.,  Vol.  11, 
pp.  894-908  (1974).  — 

7.  Jacobs,  D.A.H. : "The  Strongly  Implicit 
procedure  for  Biharmonic  Problems".  J.  of 
Comput.  Physics,  Vol.  13,  pp.  303-315 
(1973) . 


Turbulent  Flows  held  at  Univ.  College  Swan- 
sea. Proc.  to  be  published. 

16.  Taylor,  T.D.  and  Ndefo,  E.:  "Computa- 
tion of  Viscous  Flow  in  a Channel  by  the 
Method  of  Splitting".  Proc.  of  Second  In- 
ternational Conf.  on  Numerical  Methods  in 
Fluid  Dynamics,  Springer-Verlog  Publisher, 
pp.  356-364  (1970). 

17.  Kober,  H.:  "Dictionary  of  Conformal 
Representations".  Dover  Publications, 
pp.  161  (1957). 

18.  Thom,  A.:  "The  Flow  Past  Circular 
Cylinders  at  Low  Speeds".  Proc.  Roy.  Soc. 
A141,  pp.  651-666  (1933) . 

19.  Kawaguti,  M. : "Numerical  Solution  of 
the  Navier-Stokes  Equations  for  the  Flow 
Around  a Circular  Cylinder  at  Reynolds 
Number  40".  J.  Phys.  Soc.  Japan,  Vol.  8, 
pp.  747-757  (1953). 

20.  Takami,  H.  and  Keller,  H.B.:  "Steady 
Two  Dimensional  Viscous  Flow  of  an  Imcom- 
pressible  Fluid  Past  a Circular  Cylinder"/ 
Phys.  Fluids  Supplement  II,  Vol.  12, 

pp.  51-56  (1969). 

21.  Dennis,  S.C.R.  and  Chang,  G.Z.:  "Numer' 
ical  Solutions  for  Steady  Flow  Past  a 
Circular  Cylinder  at  Reynolds  Numbers  up  to 
100".  J.  Fluid  Mech.,  Vol.  42,  pp.  471- 
489  (1970) . 

22.  Lin,  C.L.,  Pepper,  D.W.  and  Lee,  S.C.: 
"Numerical  Methods  for  Separated  Flow  Solu- 
tions Around  a Circular  Cylinder".  Proc. 
AIAA  Second  Computational  Fluid  Dynamics 
Conf.  Hartford,  Connecticut,  pp.  91-100 
(1975)  . 


8.  Rubin,  S.G.  and  Khosla,  P.K.:  "Poly- 
nomial Interpolation  Methods  for  Viscous 
Plow  Calculations".  J.  Comput.  Physics. 
Vol.  24,  pp.  217-244  (1977). 

9.  Rubin,  S.G.  and  Khosla,  P.K.:  "A 
Simplified  Spline  Solution  Procedure". 

Proc.  of  6th  International  Conference  on 
Numerical  Methods  in  Fluid  Dynamics  held 
in  Tbilisi,  USSR.  (1978). 

10.  Khosla,  P.K.  and  Rubin,  S.G.:  "A  Diag- 
onally Dominant  Second-Order  Accurate  Im- 
plicit Scheme".  Computers  and  Fluids, 

Vol.  2,  pp.  207-209  (1974). 

11.  Proc.  of  1st  International  Conf.  Numer- 
ical Methods  in  Laminar  and  Turbulent 


Flows  held  at  University  College  Swansea, 
pp.  873-883  (1978). 

12.  Burggraf,  O.R. : "Analytical  and  Numer- 
ical Studies  of  Steady  Separated  Flows" . 

J.  Fluid  Mech.,  Vol.  24,  pp.  113-151  (1966). 

13.  Smith,  R.E.  and  Kl3d,  Amy:  "Compara- 


tive Study  of  Two  Numerical  Techniques  for 
the  Solution  of  Viscous  Flow  in  a Driven 
Cavity".  NASA  SP- 378  (1975). 

14.  Nallasmy,  M.  and  Krishna,  Prasad,  K.: 
"On  Cavity  Flow  at  High  Reynolds  Numbers". 
J.  Fluid  Mech.,  Vol.  2i*  PP-  391-414  (1977). 

15.  Special  Session  at  the  International 
Conf.  on  Numerical  Methods  in  Laminar  and 


G21  = " D1T1{i_l»j)-D2T5(i'l»j)  G22  = ' D3T1(i-lfj)-D4T5(i-X,j) 

G41  “ - D1T3(i-lfj)-D2T7(i-l.j)  G42  - - D3T3(i-lfj)-D4T7(i-l,j) 


G11  = " AlT2(i' j-l)-A2T6(i'j-1)  Gi2  - - A3T2(i, j-l)-A4T6(i, j-1) 

G31  * * AjT^i.j-lJ-AjTgU.j-l)  G32  - - A3T4(i,j-l)-A4T8(i,j-l) 

S11  = “i+l, j-1'  S21  = “i-l, j+1  S31  = *i+l,j-l'  S41  = *i-l,j+l 


SG1  = G1+G11S11+G21S21+G31S31+G41S41 
sg2  = g2+g12s11+g22s21+g32s31+g42s41 


GG2  = SG2-A3GN1(i, j-l)-D3GM1(i-l, j) 
-A4GM2(i, j-l)-D4GM2(i-l, j) 

Z2  = Bj+AjTjU,  j-D+DjT^i-l,  j) 

+ A2T7(i,j-l)+D2Tg(i-l, j) 

Z4  = B4+A3T3(i, j-l)+D3T4(i-l, j) 

+ A4T7(i,j-l)+D4Tg(i-l,j) 

GMX (i, j)  = (Z4GG1-Z2GG2)/DM 

tx  {i, j)  = (z2c3-z4c1)/dm 

T3(i,j)  = (Z2C4-Z4C2)/DM 
T5(i,j)  = (ZgCj-ZjCjJ/DM 
T?(i, j)  • (Z3C2-Z1C4)/DM 


GG1  = SGj^-A^Mjd,  j-l)-D1GM1(i-l,  j) 

-a2gm2 ( i , j - i ) -d2gm2 ( i- 1 , j ) 

Z1  = Bi+AiT1(i»3-1,+D1T2(i”1* j) 

+ A2T5(i, j-l)+D2T6(i-lf j) 

Z3  = Bj+AjTj  (i,  j-D+DjTj  (i-1/ j) 

+ A4T5(i,j-l)+D4T6(i-l, j) 

DM  = Z3Z4-Z2Z3 
GM2(i,j)  = (Zj^GGj-ZjGG^)  /DM 
T2(i,j)  = (Z^g-ZgE^/DN 
T4(ifj)  = (Z2E4-Z4E2)/DM 
Tg  (i,  j)  = (ZgE^-Z^Eg) /DM 
Tg  (i, j)  = (Z3E2-Z1E4)/DM 


